
      ! Simple Euler forward integrator
      ! This subroutine both advects and diffuses
      ! x0 is the initial concentration (in)
      ! x is the operated concentration (inout)
      ! u is the wind [m/s]
      ! k is the diffusion [m^2/s]
      ! n is the size of the grid 
      ! dx is the size of each grid cell [m]
      ! dt is the time step for each iteration [s]
      SUBROUTINE advf(x0, x, u, k, n, dx, dt)
         implicit none;
      
         INTEGER*4 :: i;
         INTEGER*4, INTENT(in) :: n
         INTEGER*4 :: dx, dt;
         DOUBLE PRECISION :: last_x, next_x, x_nm1, x_a, x_d;
         DOUBLE PRECISION, INTENT(in), DIMENSION(n) :: x0, u, k;
         DOUBLE PRECISION, INTENT(inout), DIMENSION(n) :: x;

         ! Loop through the grid
         do 1001 i=1,n
      
             !
             ! Advect in x
             !
             
             ! Find out what the last and next values were (considering boundary
             ! conditions)
             call get_last_val(last_x, x, i, n)
             call get_next_val(next_x, x, i, n)

             x_a = -u(i)*(next_x-last_x)/(2*dx);
      
             !
             ! Diffuse in x
             !
      
             call get_last_val(last_x, x, i, n)
             call get_next_val(next_x, x, i, n)

             x_d = k(i)*((next_x - 2*x(i) + last_x)/(dx**2));

             ! Put these together
             x(i) = dt*(x_a + x_d) + x(i);

1001  continue;

      END SUBROUTINE advf


      ! Returns the previous concentration in the array (taking boundary conditions into
      ! account)
      subroutine get_last_val(last_x, x, i, n)
          implicit none;
          INTEGER*4, INTENT(in) :: i, n;
          DOUBLE PRECISION, INTENT(inout) :: last_x;
          DOUBLE PRECISION, INTENT(in), DIMENSION(n) :: x;
          IF (i == 1) THEN
              last_x = 1.0; ! Constant influx
          ELSE IF (i == n) THEN
              last_x = x(i-1);
          ELSE
              last_x = x(i-1);
          END IF
      END SUBROUTINE get_last_val

      ! Returns the next concentration in the array (taking boundary conditions into
      ! account)
      SUBROUTINE get_next_val(next_x, x, i, n)
          implicit none;
          INTEGER*4, INTENT(in) :: i, n;
          DOUBLE PRECISION, INTENT(inout) :: next_x;
          DOUBLE PRECISION, INTENT(in), DIMENSION(n) :: x;
          IF (i == 1) THEN
              next_x = x(i+1);
          ELSE IF (i == n) THEN
              next_x = 1.0; ! Constant outflux
          ELSE
              next_x = x(i+1);
          END IF
      END SUBROUTINE get_next_val

